function [fuzzy_center, fuzzy_member]=fuzzy_nGas(x,y0,c,beta,epochs)
%introduction to parten recoginition: p/199
%for fuzzy c-means clustering it works fine for 2,4 dimension but nonsense for hingh dimension data??
%
%for neuron gas see spme competitive learning methods
%
%Input : x is the data set
%        y0 is the cluster center
%        beta is the fuzz strength
%Output: fuzzy_center
%        fuzzy_member
t0=clock;
 
Neurons=y0;
%Neurons = (rand(n,dim)-0.5)*10e-5; % small initial values
[nrow ndim]=size(y0);
D=x;
n=c; %num of neunon =num of center
%epochs=100; %num of learning need enough learning step to get good results
[dlen,dim] = size(D);
train_len = epochs*dlen;
alpha0 = 0.5; %adjust alpha_f and lmbdaf effect the results
lambda0 = 10; 

% random sample order
 rand('state',sum(100*clock));
 sample_inds = ceil(dlen*rand(train_len,1));
% lambda 
 lambda = lambda0 * (0.01/lambda0).^([0:(train_len-1)]/train_len);
% alpha
 alpha = alpha0 * (0.005/alpha0).^([0:(train_len-1)]/train_len); 
 %change alphaf from 0.005 to 0.001 effect the results

for i=1:train_len
    % sample vector
    x = D(sample_inds(i),:); % sample vector
    %x=D;
    known = ~isnan(x);       % its known components
    X = x(ones(n,1),known);  % we'll need this 
    
    
    % neighborhood ranking
    Dx = Neurons(:,known) - X;  % difference between vector and all map units
    %Dx=sqrt(som_eucdist2(Neurons(:,known),X));
    [qerrs, inds] = sort((Dx.^2)*known'); % 1-BMU, 2-BMU, etc.
    %[qerrs, inds] = sort((Dx(:,1).^2));
    ranking(inds) = [0:(n-1)];             
    h = exp(-ranking/lambda(i));
    H = h(ones(length(known),1),:)';

    % update
    oldNeurons=Neurons;
    Neurons = Neurons + alpha(i)*H.*(x(ones(n,1),known) - Neurons(:,known));
end  %end all loop

%Calculate the initial fuzzy membership based on neuron gas results
%step2 
%calculate the distance between Neuron and all other data
        d=sqrt(som_eucdist2(Neurons,D));
%step3
%calcualte the membership
        sum_c_dist=sum(d.^(-2/(beta-1)));
        membership=[d.^(-2/(beta-1))./repmat(sum_c_dist,c,1)];
%step4 updata centers
%find di=0 and replace membership as 1 all other center's membership=0 
        [zeroD_i zeroD_j]=find(Dx==0);
         membership(:,zeroD_j)=0;
        lnZero=length(zeroD_i);
        for ki=1:lnZero
            membership(zeroD_i(ki),zeroD_j(ki))=1;
        end
        fuzzy_center=Neurons;
        fuzzy_member=membership;
%step5 check error
        fuzzy_error=0;
        for ii=1:c
           fuzzy_error=fuzzy_error+som_eucdist2(Neurons(ii,:),oldNeurons(ii,:));
        end
        fuzzy_error=sqrt(fuzzy_error);
    etime(clock,t0) ;
    
    
    
    
    
    
    
